Qu’est ce que R ?
R (https://www.r-project.org/) est un langage de programmation dédié à l’analyse statistique. Libre et gratuit, disponible sur toutes les plateformes (Mac OS, Windows, Linux), il s’est largement imposé ces dernières années dans le domaine des sciences humaines et de la datascience. R dispose par défaut d’instructions de base permettant de réaliser les opérations les plus courantes. De nouvelles fonctionnalités, regroupées en packages, peuvent être facilement ajoutées pour venir enrichir les possibilités du langage. Le point fort de R par rapport à des logiciels à interfaces graphiques, est qu’il permet très facilement l’automatisation et la reproductibilité des traitements.
Environnement de travail
Rstudio (https://www.rstudio.com/) est un environnement de travail permettant d’écrire du code R. Il ne s’agit pas d’un environnement clic bouton mais d’une véritable interface de développement (IDE) destiné à écrire du code. Il s’agit en fait d’une couche logiciel permettant d’écrire, d’exécuter et de visualiser les résultats d’un programme écrit dans le langage R. RStudio n’est pas indispensable pour utiliser R. Néanmoins, il est très pratique et très bien fait.
Le menu Session > Set Working Directory permet de sélectionner le répertoire par défaut contenant vos scripts et éventuels fichiers de données. Vous pouvez également le modifier via une ligne de code :
setwd('C:/Users/Stephanie/Desktop/R_initiation')
Dans Rstudio pour compiler une ligne, vous pouvez placer le curseur dessus et cliquer sur Run, ou bien utiliser le raccourci clavier ctrl+entrée
Sous R, les éléments de base sont des objets : vecteurs, matrices, listes, table … Ces objets peuvent contenir des éléments de type numérique, booléen (logical : TRUE, FALSE), chaîne de caractère (string), facteur (factor, pour les variables qualitatives prenant un nombre déterminé de modalités). La différence entre la table de données (data.frame) et la matrice tient notamment dans le fait que la première peut contenir des éléments de types différents.
On commence par créer des éléments pour illustrer les opérations de base en R. Le signe <- permet d’assigner un contenu à une variable dont on indique le nom à gauche.
element_a <- 2
element_a
## [1] 2
element_b <- 3
element_b
## [1] 3
element_c <- element_a/2 + 2*element_b # les opérations -, *, /, sqrt pour racine, ^ pour puissance, log, exp, sont également possibles !
element_c
## [1] 7
Vecteurs
On peut travailler très naturellement sur des vecteurs
vecteur_a <- c(2,3,5)
vecteur_a
## [1] 2 3 5
On peut récupérer un élément d’un vecteur avec l’opérateur [
vecteur_a
## [1] 2 3 5
vecteur_a[1] #on commence à numéroter à partir de 1 (et non de zéro)
## [1] 2
vecteur_a[c(1,3)]
## [1] 2 5
Dans le cas où le vecteur correspond à une séquence de nombres, on peut utiliser une syntaxe particulière
vecteur_a <- c(2,3,4)
vecteur_b <- 2:4
vecteur_c <- seq(2,4,1) # les arguments de la fonction seq (séquence) correspondent à from, to, by (à partir de 2, jusqu'à 4, avec des écarts de 1)
Remarque, si on ne sait plus la signification des trois arguments de seq, on peut aller dans l’aide avec ?seq
vecteur_a == vecteur_b # compare les éléments un à un
## [1] TRUE TRUE TRUE
variable_logique_a <- vecteur_a == vecteur_b
Les opérations sur les vecteurs sont très proches des opérations sur les éléments. Une opération entre un vecteur a et un vecteur b revient à réaliser l’opération entre les couples d’éléments issus des deux vecteurs et situés à la même place (il faut que les vecteurs aient la même taille pour que ça ait un sens)
vecteur_c <- 2*vecteur_c # on écrase la valeur du vecteur_c
vecteur_c + vecteur_a/2 # exemple d'opération entre deux vecteurs
## [1] 5.0 7.5 10.0
Pour concaténer deux vecteurs, on fera tout simplement
vecteur_a
## [1] 2 3 4
vecteur_b
## [1] 2 3 4
vecteur_d <- c(vecteur_a,vecteur_b)
vecteur_d
## [1] 2 3 4 2 3 4
Le vecteur est caractérisé notamment par sa taille
length(vecteur_a)
## [1] 3
Il peut-être trié :
vecteur_d <- c(3,1,2)
order(vecteur_d)
## [1] 2 3 1
sort(vecteur_d)
## [1] 1 2 3
Listes
liste_a <- list(2,3,5)
liste_a
## [[1]]
## [1] 2
##
## [[2]]
## [1] 3
##
## [[3]]
## [1] 5
liste_a[[2]]
## [1] 3
liste_b <- list(sexe = 2, age = 3, salaire = 5)
liste_b$salaire # si on donne des noms aux éléments de la liste, on peut les récupérer via l'opérateur $
## [1] 5
Pour concaténer deux listes
liste_a
## [[1]]
## [1] 2
##
## [[2]]
## [1] 3
##
## [[3]]
## [1] 5
liste_a <- append(liste_a, element_a) # ou liste_b à la place de element_a pour concaténer deux listes
liste_a
## [[1]]
## [1] 2
##
## [[2]]
## [1] 3
##
## [[3]]
## [1] 5
##
## [[4]]
## [1] 2
La fonction length convient aussi pour les listes
length(liste_a)
## [1] 4
on peut passer de la liste au vecteur avec la fonction unlist
unlist(liste_a)
## [1] 2 3 5 2
Matrices
matrice_a <- matrix(1:15,ncol=5) # exemple d'une matrice remplie des chiffres consécutifs de 1 à 15 rangés sur 5 colonnes
matrice_a
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 4 7 10 13
## [2,] 2 5 8 11 14
## [3,] 3 6 9 12 15
head(matrice_a) # la fonction head permet de visualiser un extrait des données, ici elles sont petites donc c'est la totalité
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 4 7 10 13
## [2,] 2 5 8 11 14
## [3,] 3 6 9 12 15
matrice_a[1,2] # pour récupérer l'élément de la première ligne et de la deuxième colonne
## [1] 4
Par défaut, les opérations mathématiques simples se font élément par élément
matrice_b <- 3*matrice_a
matrice_a
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 4 7 10 13
## [2,] 2 5 8 11 14
## [3,] 3 6 9 12 15
matrice_b
## [,1] [,2] [,3] [,4] [,5]
## [1,] 3 12 21 30 39
## [2,] 6 15 24 33 42
## [3,] 9 18 27 36 45
matrice_b - matrice_a # les opérations se font éléments par élément
## [,1] [,2] [,3] [,4] [,5]
## [1,] 2 8 14 20 26
## [2,] 4 10 16 22 28
## [3,] 6 12 18 24 30
t(matrice_b)%*% matrice_a # pour faire des vrais produits matriciels, on utilise l'opérateur %*%, ici t() indique que l'on prend également la transposée
## [,1] [,2] [,3] [,4] [,5]
## [1,] 42 96 150 204 258
## [2,] 96 231 366 501 636
## [3,] 150 366 582 798 1014
## [4,] 204 501 798 1095 1392
## [5,] 258 636 1014 1392 1770
Les dimensions de la matrice peuvent être obtenues de la manière suivante :
nrow(matrice_a)
## [1] 3
ncol(matrice_a)
## [1] 5
dim(matrice_a)
## [1] 3 5
Pour concaténer deux matrices, on peut soit les mettre côte à côte (cbind), soit les empiler (rbind). Ces fonctions fonctionnent aussi pour les data.frames.
matrice_c <- cbind(matrice_a, matrice_b)
matrice_c
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1 4 7 10 13 3 12 21 30 39
## [2,] 2 5 8 11 14 6 15 24 33 42
## [3,] 3 6 9 12 15 9 18 27 36 45
dim(matrice_c)
## [1] 3 10
matrice_d <- rbind(matrice_a, matrice_b)
matrice_d
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 4 7 10 13
## [2,] 2 5 8 11 14
## [3,] 3 6 9 12 15
## [4,] 3 12 21 30 39
## [5,] 6 15 24 33 42
## [6,] 9 18 27 36 45
dim(matrice_d)
## [1] 6 5
Data.frame
Les tables de données ou data.frame est sans doute le format qu’on utilisera le plus dans la suite.
df_a <- as.data.frame(matrice_a) # on commence par transformer la matrice en data.frame, l'opération symétrique as.matrix() est possible
head(df_a)
names(df_a) # ici le data.frame provient de la conversion d'une matrice, il n'y a donc pas de noms de colonnes à part ceux qui sont donnés par défaut V1, V2...
## [1] "V1" "V2" "V3" "V4" "V5"
names(df_a) <- c('a','b','c','d','e') # on peut changer le nom des colonnes simplement
Pour récupérer de l’information dans la table de données, on peut procéder de différentes manières
df_a$b # avec l'opérateur $, on peut récupérer directement la colonne b
## [1] 4 5 6
df_a[['b']] # revient au même
## [1] 4 5 6
df_a[[2]] # correspond également car il s'agit de la 2e colonne
## [1] 4 5 6
df_a$e <- 1 # permet de construire une nouvelle colonne e remplie de 1
head(df_a)
df_a$f <- df_a$a + df_a$b # permet de construire une nouvelle colonne qui serait la somme des deux premières, les opération -, *, /, sqrt pour racine carrée, log, exp, ^, sont également possibles !
head(df_a)
df_a$f <- sqrt(df_a$a) + df_a$b^2 # permet de construire une nouvelle colonne qui serait la somme des deux premières, les opération -, *, /, sqrt pour racine
Vérifions le type de ces objets avec la fonction class
class(element_a)
## [1] "numeric"
class(vecteur_a)
## [1] "numeric"
class(variable_logique_a)
## [1] "logical"
class(liste_a)
## [1] "list"
class(matrice_a)
## [1] "matrix"
class(df_a)
## [1] "data.frame"
Nous avons créé un certain nombre de variables qui sont disponibles dans votre environnement (généralement en haut à droite par défaut) On peut aussi utiliser la fonction ls pour voir la liste des objets dont on dispose, et rm pour en supprimer.
ls()
## [1] "df_a" "element_a" "element_b"
## [4] "element_c" "liste_a" "liste_b"
## [7] "matrice_a" "matrice_b" "matrice_c"
## [10] "matrice_d" "variable_logique_a" "vecteur_a"
## [13] "vecteur_b" "vecteur_c" "vecteur_d"
rm(element_a)
ls()
## [1] "df_a" "element_b" "element_c"
## [4] "liste_a" "liste_b" "matrice_a"
## [7] "matrice_b" "matrice_c" "matrice_d"
## [10] "variable_logique_a" "vecteur_a" "vecteur_b"
## [13] "vecteur_c" "vecteur_d"
Reprenons le cas où on a un data.frame simple. On veut calculer la moyenne, la somme, le maximum et le minimum de chaque colonne.
Une manière de le faire serait de faire une boucle sur le nombre de colonnes et d’appliquer la fonction mean, sum, max, min, à chaque colonne tour à tour.
nb_col <- ncol(df_a)
mean_col <- NULL # on commence par créer une variable vide, dans laquelle on va ajouter itérativement les moyenne de chaque colonne
for (i in 1:nb_col){
print(i)
mean_temp <- mean(df_a[[i]])
print(mean_temp)
mean_col <- c(mean_col, mean_temp)
print(mean_col)
print('---------')
}
## [1] 1
## [1] 2
## [1] 2
## [1] "---------"
## [1] 2
## [1] 5
## [1] 2 5
## [1] "---------"
## [1] 3
## [1] 8
## [1] 2 5 8
## [1] "---------"
## [1] 4
## [1] 11
## [1] 2 5 8 11
## [1] "---------"
## [1] 5
## [1] 1
## [1] 2 5 8 11 1
## [1] "---------"
## [1] 6
## [1] 27.04875
## [1] 2.00000 5.00000 8.00000 11.00000 1.00000 27.04875
## [1] "---------"
En réalité, les boucles sont à éviter en R, elles ne sont pas efficaces, on leur préférera la fonction lapply qui distribue la fonction souhaitée sur chaque élément d’une liste (et retourne une liste). Un data.frame peut-être vu comme une liste de vecteurs.
mean_col <- lapply(df_a, mean) # ici la fonction est particulièrement simple car mean est déjà une fonction R
mean_col <- sapply(df_a, mean) # sapply procède de même mais renvoie un vecteur
mean_col
## a b c d e f
## 2.00000 5.00000 8.00000 11.00000 1.00000 27.04875
Admettons qu’on cherche à appliquer une fonction personnalisée
mean_personnalisee <- function(vect){
return(sum(vect)/length(vect))
}
mean_col <- sapply(df_a, function(vect) mean_personnalisee(vect))
mean_col
## a b c d e f
## 2.00000 5.00000 8.00000 11.00000 1.00000 27.04875
Appliquée sur un vecteur, lapply le convertit en liste, on peut repasser à un format vecteur avec la fonction unlist
fonction_personnalisee <- function(x){
return(x*log(x))
}
vecteur_a_transforme <- lapply(vecteur_a, function(x) fonction_personnalisee(x))
class(vecteur_a_transforme)
## [1] "list"
vecteur_a_transforme
## [[1]]
## [1] 1.386294
##
## [[2]]
## [1] 3.295837
##
## [[3]]
## [1] 5.545177
unlist(vecteur_a_transforme)
## [1] 1.386294 3.295837 5.545177
vecteur_a*log(vecteur_a) # aurait marché aussi !
## [1] 1.386294 3.295837 5.545177
apply s’applique aussi à des matrices, on indique alors si la fonction à distribuer doit être distribuée en ligne ou en colonne
matrice_a
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 4 7 10 13
## [2,] 2 5 8 11 14
## [3,] 3 6 9 12 15
apply(matrice_a,1,sum) # ici on applique la fonction somme, autrement dit on somme les éléments par ligne
## [1] 35 40 45
apply(matrice_a,2,sum) # par colonne
## [1] 6 15 24 33 42
La syntaxe de la condition est très proche de celle de la boucle, sauf que l’on remplace for par if
fonction_personnalisee <- function(x){
if (x>0){res <- x*log(x)}
else {res <- 0}
return(res)
}
vecteur_d <- c(0,-1,2)
sapply(vecteur_d, function(x) fonction_personnalisee(x))
## [1] 0.000000 0.000000 1.386294
abs(vecteur_d) # rmq : abs permet de passer un élément en valeur absolue (s'applique également aux vecteurs sans recours à lapply)
## [1] 0 1 2
Ici on a inclus la condition dans une fonction, mais ce n’est pas indispensable ! On peut vouloir aussi tester la différence != plutôt que l’égalité, < ou >, et si on a plusieurs conditions, on mettra chacune entre parenthèse et on utilisera | pour dire “ou” et & pour dire “et” (par exemple (b-3==a) & (b>=a)).
Parfois, on n’a pas besoin de passer par if pour appliquer une condition, par exemple si on veut récupérer seulement les éléments de vecteur_d supérieurs ou égaux à 0. Cela revient à dire qu’on teste la condition sur chaque élément et qu’on conserve vecteur_d[i] pour i tel que vecteur_d[i]>=0
vecteur_d
## [1] 0 -1 2
vecteur_d[vecteur_d>=0]
## [1] 0 2
Et ça fonctionne aussi avec les data.frames
df_a[df_a$a>2] # on ne garde que les colonnes qui vérifient cette condition
Il arrive que les objets que l’on manipule ne soit pas numériques, on parle alors de character ou de factor. C’est le cas par exemple des noms de colonnes de df_a
class(names(df_a))
## [1] "character"
On peut aussi réaliser des opérations, par exemple de concaténation. Imaginons que l’on souhaite préciser les noms de colonnes en indiquant qu’il s’agit de département.
paste('departement','test',sep='_')
## [1] "departement_test"
paste(rep('departement',5), names(df_a),sep='_')
## [1] "departement_a" "departement_b" "departement_c" "departement_d"
## [5] "departement_e" "departement_f"
names(df_a) <- paste(rep('departement',5), names(df_a),sep='_') # revient à faire apply de paste sur tous les éléments du vecteur names(df_a)
rep('departement',5) # rep est la fonction qui permet de créer des vecteurs d'une taille donnée contenant toujours le même élément qui pourrait être numérique
## [1] "departement" "departement" "departement" "departement" "departement"
On peut récupérer une partie d’une chaine de caractère simplement
character_a <- 'departement_a'
class(character_a)
## [1] "character"
substr(character_a,1,10) #les arguments correspondent à la chaine, au début de la sous-chaîne que l'on veut récupérer, et à la fin de la sous-chaîne que l'on veut récupérer
## [1] "departemen"
strsplit(character_a,'_') #les arguments correspondent à la chaine et au caractère qui au niveau duquel on veut couper
## [[1]]
## [1] "departement" "a"
strsplit(character_a,'_')[[1]][[2]] #le résultat est une liste donc on utilise les doubles crochets pour récupérer l'élément qui nous intéresse
## [1] "a"
La taille d’une chaine de caractères correspond au nombre de caractères
nchar(character_a)
## [1] 13
Les facteurs correspondent aux modalités d’une variable qualitative.
col_df_a <- names(df_a)
col_df_a_fact <- as.factor(col_df_a)
is.character(col_df_a_fact)
## [1] FALSE
class(col_df_a)
## [1] "character"
class(col_df_a_fact)
## [1] "factor"
On ne peut pas rajouter un élément qui n’est pas dans la liste des facteurs naïvement
levels(col_df_a_fact)
## [1] "departement_a" "departement_b" "departement_c" "departement_d"
## [5] "departement_e" "departement_f"
col_df_a_fact <- c(col_df_a_fact, 'departement_z')
class(col_df_a_fact) # attention cette opération a converti col_df_a_fact en character
## [1] "character"
col_df_a_fact
## [1] "1" "2" "3" "4"
## [5] "5" "6" "departement_z"
levels(col_df_a_fact)
## NULL
col_df_a_fact <- as.factor(col_df_a)
as.character(col_df_a_fact) # parfois on préfère revenir aux chaines de caractères pour éviter ce type de problèmes
## [1] "departement_a" "departement_b" "departement_c" "departement_d"
## [5] "departement_e" "departement_f"
Lorsque l’on charge des tables de données, il est possible que R ne reconnaisse pas une variable numérique et qu’elle soit codée en facteurs. Dans ce cas, il faut la convertir en chaines de caractères avant de la convertir en nombre, car les facteurs sont associés à des numéros qui correspondent aux modalités.
nombres <- c(4,2,3)
facteurs <- as.factor(nombres)
facteurs
## [1] 4 2 3
## Levels: 2 3 4
as.numeric(facteurs)
## [1] 3 1 2
as.numeric(as.character(facteurs))
## [1] 4 2 3
On a souvent besoin de générer des nombres aléatoires, par exemple pour tirer un échantillon dans les lignes d’une table, le plus simple est d’utiliser sample :
sample(1:100,5) # on tire 5 éléments sans remise entre 1 et 100, l'argument replace = TRUE permet de faire un tirage avec remise.
## [1] 4 78 27 83 97
Il est parfois intéressant de simuler des données. Générons 50 observations tirées d’une variable suivant une loi normale de moyenne 20 et d’écart-type 10.
set.seed(123) # permet de reproduire les mêmes résultats d'une fois sur l'autre en dépit de l'aléa présent
n <- 50
Y <- rnorm(n,20,5)
mean(Y) # moyenne empirique
## [1] 20.17202
sd(Y) # écart-type empirique
## [1] 4.62935
sd(Y)/sqrt(length(Y)) # écart-type de la moyenne empirique
## [1] 0.6546889
summary(Y) # quartiles et moyenne empiriques
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.17 17.20 19.64 20.17 23.49 30.84
On peut représenter graphiquement cette variable
boxplot(Y) # diagramme boîte
hist(Y, probability=T, col="blue") # histogramme de la densité
lines(density(Y), col="red", lwd=2) # lissage de l'histogramme
# tracer la loi théorique
x <- 1:100
curve(dnorm(x,mean=20,sd=5),add=TRUE,col="green",lwd=2) # l'argument add permet de rajouter cette courbe sur le même graphique, il s'agit de la distribution théorique
On peut générer des observations tirées dans une loi uniforme sur [a,b]
X <- runif(n,50,100)
hist(X)
Vous pouvez augmenter la valeur de n pour voir ce que ça donne
Pour cette entrée en matière de l’exploration d’un vrai fichier de donnée, nous nous sommes fortement inspirés de l’exercice suivant https://www.math.univ-toulouse.fr/~besse/Wikistat/pdf/st-scenar-statlab.pdf
Une étude réalisée entre 1961 et 1973 dans la maternité d’un hôpital d’Oakland (Californie) avait pour but de rechercher si certaines caractéristiques des parents avaient une influence sur le développement de l’enfant. Parmi les variables collectées, 19 variables décrites dans le tableau ci-dessous ont été observées sur 115 familles ou unités statistiques. Ces variables décrivent des informations médicales et socio-économiques concernant le bébé et ses parents au moment de la naissance puis dix ans plus tard, permettant ainsi de se poser différentes questions de nature plutôt épidémiologique :
famil=read.csv2("statlab.csv") # lecture du fichier csv
head(famil) # aperçu du haut des données
dim(famil) # combien de colonnes et de ligne
## [1] 115 19
summary(famil) # quelques statistiques classiques
## ESx ERh ET0 EP0 ET10
## F:46 R-:22 Min. :43.00 Min. :1.770 Min. :124.0
## M:69 R+:93 1st Qu.:51.00 1st Qu.:2.925 1st Qu.:132.5
## Median :52.00 Median :3.360 Median :136.0
## Mean :51.93 Mean :3.389 Mean :136.0
## 3rd Qu.:53.50 3rd Qu.:3.760 3rd Qu.:139.8
## Max. :58.50 Max. :6.350 Max. :151.5
## EP10 MRh MA0 MP0 MCig0
## Min. :23.10 R-:19 Min. :15.00 Min. : 43.50 >10cig :39
## 1st Qu.:27.45 R+:96 1st Qu.:24.00 1st Qu.: 51.50 0cig :48
## Median :31.80 Median :27.00 Median : 56.50 1-10cig:28
## Mean :32.53 Mean :27.53 Mean : 59.57
## 3rd Qu.:36.05 3rd Qu.:32.00 3rd Qu.: 63.25
## Max. :57.20 Max. :43.00 Max. :113.50
## MT MP10 MCig10 PA0
## Min. :148.5 Min. : 38.50 >10cig :39 Min. :20.00
## 1st Qu.:160.2 1st Qu.: 55.25 0cig :55 1st Qu.:26.00
## Median :164.0 Median : 61.00 1-10cig:21 Median :29.00
## Mean :164.1 Mean : 64.07 Mean :30.57
## 3rd Qu.:167.8 3rd Qu.: 69.25 3rd Qu.:35.00
## Max. :178.0 Max. :115.50 Max. :50.00
## PCig0 PT PP10 RF0
## >10cig :65 Min. :158.0 Min. : 49.50 Min. : 14.0
## 0cig :36 1st Qu.:172.8 1st Qu.: 72.50 1st Qu.: 54.0
## 1-10cig:14 Median :179.0 Median : 80.50 Median : 72.0
## Mean :179.0 Mean : 82.23 Mean : 76.1
## 3rd Qu.:185.5 3rd Qu.: 90.50 3rd Qu.: 93.0
## Max. :198.0 Max. :129.00 Max. :250.0
## RF10
## Min. : 21.0
## 1st Qu.:110.0
## Median :150.0
## Mean :156.5
## 3rd Qu.:190.0
## Max. :500.0
La fonction summary fournit déjà beaucoup d’information sur chaque variable. Mais on pourrait réappliquer les fonctions moyenne, écart-type, quantiles, diagramme boîte, histogramme, indépendamment, sur toutes les variables ou certaines en particulier.
sapply(famil, mean) # les moyennes de chaque variable
## ESx ERh ET0 EP0 ET10 EP10
## NA NA 51.930435 3.388783 136.047826 32.526087
## MRh MA0 MP0 MCig0 MT MP10
## NA 27.530435 59.573913 NA 164.104348 64.065217
## MCig10 PA0 PCig0 PT PP10 RF0
## NA 30.565217 NA 178.965217 82.234783 76.095652
## RF10
## 156.495652
sapply(famil, function(x) mean(x,na.rm = T)) # les moyennes de chaque variable en ignorant les valeurs manquantes
## ESx ERh ET0 EP0 ET10 EP10
## NA NA 51.930435 3.388783 136.047826 32.526087
## MRh MA0 MP0 MCig0 MT MP10
## NA 27.530435 59.573913 NA 164.104348 64.065217
## MCig10 PA0 PCig0 PT PP10 RF0
## NA 30.565217 NA 178.965217 82.234783 76.095652
## RF10
## 156.495652
sapply(famil, sd) # les écarts type
## ESx ERh ET0 EP0 ET10 EP10
## 0.4920419 0.3950495 2.9793787 0.6767183 6.0244290 6.1476438
## MRh MA0 MP0 MCig0 MT MP10
## 0.3730019 5.8524575 12.9283185 0.7605851 5.5025894 13.3862200
## MCig10 PA0 PCig0 PT PP10 RF0
## 0.7082385 6.4673681 0.7032669 7.7521140 13.6865476 36.1429845
## RF10
## 72.5150165
boxplot(famil$ET0) # taille de l'enfant
boxplot(famil$EP0) # poids de l'enfant
hist(famil$ET0)
hist(famil$EP0)
Pour les variables qualitatives, on aura plutôt tendance à considérer les fréquences des modalités.
table(famil$ESx)
##
## F M
## 46 69
barplot(table(famil$ESx)) # sexe de l'enfant
barplot(table(famil$MCig0)) # consommation de cigarettes
pie(table(famil$MCig10)) # consommation de cigarettes dix ans après
Pour analyser grossièrement les liaisons, on peut regarder un nuage de points
pairs(famil[,c(3:6,8,9,11)])
plot(EP10~PP10,data=famil) # poids de l'enfant à 10 ans, poids du père
plot(EP10~ET10,data=famil) # poids de l'enfant à 10 ans, taille à 10 ans
Pour analyser grossièrement les liaisons entre variables qualitatives, on aura recours aux tables de contingence
table(famil$ESx,famil$ERh) # sexe et rhésus
##
## R- R+
## F 9 37
## M 13 56
addmargins(table(famil$ESx,famil$ERh)) # avec les marges
##
## R- R+ Sum
## F 9 37 46
## M 13 56 69
## Sum 22 93 115
prop.table(table(famil$ESx,famil$ERh)) # fréquences relatives
##
## R- R+
## F 0.07826087 0.32173913
## M 0.11304348 0.48695652
Une manière de représenter graphiquement ces tables est d’utiliser une moisaique
mosaicplot(table(famil$ESx,famil$ERh))
addmargins(table(famil$MCig0,famil$ESx)) # consommation de cigarette et sexe de l'enfant
##
## F M Sum
## >10cig 11 28 39
## 0cig 27 21 48
## 1-10cig 8 20 28
## Sum 46 69 115
mosaicplot(table(famil$MCig0,famil$ESx))
Si l’on veut croiser variables qualitatives et quantitatives, on peut utiliser les boites
boxplot(EP0~ESx,data=famil) # poids de l'enfant vs sexe
boxplot(EP0~MCig0,data=famil) # poids de l'enfant vs consommation de cigarettes
Remarque, les modalités de Mcig0 ne sont pas dans l’ordre le plus naturel, on peut décider de réordonner les modalités.
famil$MCig0 <- factor(famil$MCig0, levels = c("0cig", "1-10cig", ">10cig"))
famil$MCig10 <- factor(famil$MCig10, levels = c("0cig", "1-10cig", ">10cig"))
boxplot(EP0~MCig0,data=famil) # poids de l'enfant vs consommation de cigarettes
Quelques tests, imaginons qu’on veuille tester l’indépendance de deux variables qualitatives. Le test du \(\chi^2\) sera adapté à ce problème (dans le cas où les effectifs d’une modalité sont trop faibles, il faut regrouper les modalités).
chisq.test(table(famil$ESx,famil$ERh)) # Sexe de l'enfant vs Rhésus
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(famil$ESx, famil$ERh)
## X-squared = 5.6549e-31, df = 1, p-value = 1
chisq.test(table(famil$ESx,famil$MCig0)) # Sexe de l'enfant vs consommation de cigarettes
##
## Pearson's Chi-squared test
##
## data: table(famil$ESx, famil$MCig0)
## X-squared = 9.0657, df = 2, p-value = 0.01075
On s’intéresse à l’influence du sexe sur la taille à la naissance. Tester l’égalité des deux moyennes nécessite de vérifier préalablement plusieurs points : la normalité des distributions dans chaque classe à moins que l’échantillon soit considéré de taille suffisamment grande, le caractère indépendant ou appariés des échantillons, l’égalité ou non des variances à l’intérieur de chaque groupe. On dispose de deux échantillons indépendants : les garçons et les filles. Testons les autres hypothèses.
# Normalité des distributions (facultatif car n grand ici)
shapiro.test(famil[famil$ESx=="M","ET0"])
##
## Shapiro-Wilk normality test
##
## data: famil[famil$ESx == "M", "ET0"]
## W = 0.95676, p-value = 0.01779
shapiro.test(famil[famil$ESx=="F","ET0"])
##
## Shapiro-Wilk normality test
##
## data: famil[famil$ESx == "F", "ET0"]
## W = 0.95141, p-value = 0.05315
# égalité des variances (test de Fisher)
var.test(ET0~ESx,data=famil)
##
## F test to compare two variances
##
## data: ET0 by ESx
## F = 0.92761, num df = 45, denom df = 68, p-value = 0.7979
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.5495201 1.6132113
## sample estimates:
## ratio of variances
## 0.9276114
Le test de comparaison des moyennes à utiliser (Student vs. Welsh) dépend du résultat précédent concernant l’égalité des variances.
t.test(ET0~ESx,var.equal=F, data=famil) # si les variances sont différentes c'est un test de Welsh
##
## Welch Two Sample t-test
##
## data: ET0 by ESx
## t = 2.9017, df = 99.064, p-value = 0.004573
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.5006535 2.6660131
## sample estimates:
## mean in group F mean in group M
## 52.88043 51.29710
t.test(ET0~ESx,var.equal=T, data=famil) # Dans le cas où elles sont considérées égales, c'est un test de Student.
##
## Two Sample t-test
##
## data: ET0 by ESx
## t = 2.8798, df = 113, p-value = 0.004761
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.4940799 2.6725868
## sample estimates:
## mean in group F mean in group M
## 52.88043 51.29710
Dans le cas d’échantillons appariés, par exemple si on se propose d’étudier l’évolution du poids de la mère au moment de la naissance et dix ans après, on utilise l’option paired=TRUE
t.test(famil$MP0, famil$MP10,paired=TRUE)
##
## Paired t-test
##
## data: famil$MP0 and famil$MP10
## t = -7.458, df = 114, p-value = 1.864e-11
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -5.684285 -3.298324
## sample estimates:
## mean of the differences
## -4.491304
Remarques :
si l’hypothèse de normalité des distributions n’est pas vérifiée et si l’échantillon est trop réduit, c’est un test non-paramétrique qu’il faut mettre en ouvre. Les tests non-paramétriques sont basés sur les rangs des observations et donc sur les comparaisons des médianes des échantillons (wilcoxson).
Si l’on veut tester l’indépendance entre une variable qualitative et quantitative, l’ANOVA associée à un test de Fisher est sans doute le test le plus utilisé ; il revient au test de Student lorsque la variable qualitative n’a que deux modalités.
La régression simple permet de tester l’influence éventuelle d’une variable sur une autre et plus précisément, dans le cas de cet exemple, d’expliquer la taille de l’enfant à 10 ans. On estime le modèle avec la fonction lm
res1.reg <- lm(ET10 ~ PT, data = famil)
names(res1.reg) # liste des résultats
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"
Les graphiques suivants permettent de s’assurer de la validité du modèle, et notamment de statuer sur l’homoscédasticité des résidus, leur normalité, la bonne linéarité du modèle.
qqnorm(res1.reg$residuals) # normalité des résidus
qqline(res1.reg$residuals)
shapiro.test(res1.reg$residuals)
##
## Shapiro-Wilk normality test
##
## data: res1.reg$residuals
## W = 0.98511, p-value = 0.2346
Les résidus sont “grands” si, une fois normalisés ou plutôt “studentisés”, ils sont de valeur absolue plus grande que 2.
res.student <- rstudent(res1.reg)
ychap=res1.reg$fitted.values
plot(res.student~ychap,ylab="Résidus")
abline(h=c(-2,0,2),lty=c(2,1,2)) # rajoute une ligne horizontale
Une observation est influente si elle a un grand résidu est est associée à une grande valeur sur la diagonale de la hat matrix. Cela correspond à une valeur élevée (plus grande que 1) de la distance de Cook.
cook <- cooks.distance(res1.reg) # repérage de points influents
plot(cook~ychap,ylab="Distance de Cook")
abline(h=c(0,1),lty=c(1,2))
La significativité du modèle peut être analysée via la fonction summary :
summary(res1.reg)
##
## Call:
## lm(formula = ET10 ~ PT, data = famil)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.2297 -3.8156 -0.1753 3.4956 12.5815
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 92.82881 12.44735 7.458 1.94e-11 ***
## PT 0.24149 0.06949 3.475 0.000725 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.751 on 113 degrees of freedom
## Multiple R-squared: 0.09657, Adjusted R-squared: 0.08857
## F-statistic: 12.08 on 1 and 113 DF, p-value: 0.0007246
La régression linéaire simple conduit à un modèle très mal ajusté. Le modèle linéaire multiple ci-dessous, plus complexe, recherche un meilleur ajustement des données.
res2.reg <- lm(ET10~ET0+EP0+MA0+MP0+MT+MP10+PA0+PT+PP10+RF0+RF10, data = famil)
plot(res2.reg) # diagnostics des résidus
summary(res2.reg) # résultats
##
## Call:
## lm(formula = ET10 ~ ET0 + EP0 + MA0 + MP0 + MT + MP10 + PA0 +
## PT + PP10 + RF0 + RF10, data = famil)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.753 -3.605 -0.583 3.575 13.055
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.395052 22.732841 0.853 0.39554
## ET0 0.603298 0.293976 2.052 0.04269 *
## EP0 -1.024211 1.330923 -0.770 0.44333
## MA0 -0.070042 0.149033 -0.470 0.63937
## MP0 -0.066802 0.082435 -0.810 0.41961
## MT 0.346222 0.107974 3.207 0.00179 **
## MP10 0.078521 0.079884 0.983 0.32794
## PA0 0.147236 0.130401 1.129 0.26148
## PT 0.154598 0.089309 1.731 0.08644 .
## PP10 0.035471 0.050199 0.707 0.48140
## RF0 -0.007751 0.016938 -0.458 0.64819
## RF10 -0.010480 0.008614 -1.217 0.22652
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.245 on 103 degrees of freedom
## Multiple R-squared: 0.3152, Adjusted R-squared: 0.242
## F-statistic: 4.309 on 11 and 103 DF, p-value: 2.792e-05
On peut aller plus loin dans l’exploration sur les données sans fixer une variable cible, en réalisant une analyse en composantes principales. On ne garde que les variables quantitatives.
Pour cela, on va faire appel à la librairie prcomp, ce n’est pas la seule qui permet de faire une ACP. Assez régulièrement sous R, on va avoir besoin de mobiliser des packages ou librairies. La liste complète des packages disponibles gratuitement est consultable sur le site du CRAN. L’installation d’un package supplémentaire peut se faire via le menu Packages>Installer le(s) package(s) et en choisissant un site miroir du CRAN. On peut également télécharger l’archive .zip correspondant au package et utiliser ensuite Packages/Installer depuis des fichiers zip
On peut sinon taper la commande :
#install.packages('prcomp') # décommenter pour installer
data <- famil[,c(3:6,8,9,11,12,14,16:19)] # liste des varaibles quantitatives
noms <- names(data)
res.pca <- prcomp(data,scale=T)
plot(res.pca) # décroissance des valeurs propres
summary(res.pca) # parts de variance expliquée
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 1.8915 1.4792 1.3594 1.13565 1.04918 0.95572
## Proportion of Variance 0.2752 0.1683 0.1422 0.09921 0.08468 0.07026
## Cumulative Proportion 0.2752 0.4435 0.5857 0.68488 0.76956 0.83982
## PC7 PC8 PC9 PC10 PC11 PC12
## Standard deviation 0.82493 0.66919 0.56744 0.43219 0.4281 0.39240
## Proportion of Variance 0.05235 0.03445 0.02477 0.01437 0.0141 0.01184
## Cumulative Proportion 0.89217 0.92661 0.95138 0.96575 0.9798 0.99169
## PC13
## Standard deviation 0.32870
## Proportion of Variance 0.00831
## Cumulative Proportion 1.00000
plot(res.pca$x,col=c('black','red','green')[famil$MCig0]) # les observations sont représentées dans le plan des deux composantes principales résumant le mieux l'information contenue dans les données, la couleur permet de distinguer la consommation de cigarette (noir pour la consommation la plus élevée)
text(10*res.pca$rotation,noms,col="blue") # on peut afficher également les variables pour interpréter les axes, on voit que l'axe vertical rend surtout compte des revenus tandis que l'axe horizontal est plus lié à la taille et au poids de l'enfant et des parents.
abline(h=0,v=0,lty=2)
Cet exercice mobilise l’enquête nationale sur les structures des urgences hospitalières 2013. Un jour donné, le 11 juin 2013 de 8h à 8h le lendemain, un questionnaire papié renseigné en temps réel en parallèle de la prise en charge par tous les points d’accueils (736) des établissements de santé autorisés pour l’activité d’accueil et de traitement des urgences (634) et pour tous les patients (52018).
urgences <- read.csv2("enq_urgences_structure.csv", sep=',') # lecture du fichier csv
head(urgences) # aperçu du début des données
dim(urgences) # combien de colonnes et de lignes
## [1] 734 223
#names(urgences)
dim(urgences)
## [1] 734 223
Remarque, si le fichier était au format xls, on pourrait l’ouvrir, par exemple, avec le package xlsx
library("xlsx")
urgences2 <- read.xlsx("enq_urgences_structure.xls",1) # 1 correspond ici à la page 1 de l'excel
head(urgences2) # aperçu du haut des données
dim(urgences2)
## [1] 734 225
Pour chaque heure de pointage (8h, 12h, 18h, 22h, 8h) indiquée par les lettres A, B, C, D, E, on retient pour l’exercice les variables suivantes (cf questionnaire dans le dossier de la formation) :
Les variables de présence de patients aux urgences :
Les variables d’organisation du travail :
Les variables liées à l’établissement :
A103 : Nombre de passages sur les 24 heures A4 : type d’accueil A2 : numéro finess de l’établissement
Pour ne garder que ces variables là, on reconstruit les noms des variables pour chaque heure de pointage. Par exemple, la variable 98 : nombre d’aide soignant, va être présente cinq fois : “ID_A87”, “ID_B87”, “ID_C87”, “ID_D87”, “ID_E87”. On va utiliser lapply et paste (qui concatène les chaines de caractères) pour reconstituer tous les noms de variables qu’on veut garder pour toutes les heures de pointages. rep est une fonction qui répète plusieurs fois le même élément.
col_a_garder <- c(87:90,92:101)
nbcol <- length(col_a_garder)
paste(rep('ID_',nbcol),rep('A',nbcol),col_a_garder,sep='') # génère la liste des variables pour la première heure de pointage
## [1] "ID_A87" "ID_A88" "ID_A89" "ID_A90" "ID_A92" "ID_A93" "ID_A94"
## [8] "ID_A95" "ID_A96" "ID_A97" "ID_A98" "ID_A99" "ID_A100" "ID_A101"
col_a_garder <- unlist(lapply(c('A','B','C','D','E'), function(lettre) paste(rep('ID_',nbcol),rep(lettre,nbcol),col_a_garder,sep=''))) # génère les 5 listes de variables en une seule ligne de code
col_a_garder
## [1] "ID_A87" "ID_A88" "ID_A89" "ID_A90" "ID_A92" "ID_A93" "ID_A94"
## [8] "ID_A95" "ID_A96" "ID_A97" "ID_A98" "ID_A99" "ID_A100" "ID_A101"
## [15] "ID_B87" "ID_B88" "ID_B89" "ID_B90" "ID_B92" "ID_B93" "ID_B94"
## [22] "ID_B95" "ID_B96" "ID_B97" "ID_B98" "ID_B99" "ID_B100" "ID_B101"
## [29] "ID_C87" "ID_C88" "ID_C89" "ID_C90" "ID_C92" "ID_C93" "ID_C94"
## [36] "ID_C95" "ID_C96" "ID_C97" "ID_C98" "ID_C99" "ID_C100" "ID_C101"
## [43] "ID_D87" "ID_D88" "ID_D89" "ID_D90" "ID_D92" "ID_D93" "ID_D94"
## [50] "ID_D95" "ID_D96" "ID_D97" "ID_D98" "ID_D99" "ID_D100" "ID_D101"
## [57] "ID_E87" "ID_E88" "ID_E89" "ID_E90" "ID_E92" "ID_E93" "ID_E94"
## [64] "ID_E95" "ID_E96" "ID_E97" "ID_E98" "ID_E99" "ID_E100" "ID_E101"
col_a_garder <- c(col_a_garder, "ID_A103","ID_A4","ID_A2")
urgences <- urgences[,col_a_garder]
dim(urgences)
## [1] 734 73
Regardons la répartition des types d’établissement
table(urgences$ID_A4) # répartition des structures
##
## Structures des urgences générales
## 536
## Structures des urgences générales sur site ayant aussi une structure des urgences pédiatriques
## 94
## Structures des urgences pédiatriques
## 104
On va regarder la distribution du nombre de passages total dans les établissements, c’est un proxy de la taille de ces derniers. On va créer des modalités pour identifier les services d’urgences comparables. Pour les graphiques, on va désormais utiliser la librairie ggplot2, installez-là si elle n’est pas installée. aes permet de définir les variables qui vont être utilisées.
library(ggplot2)
ggplot(urgences, aes(x=ID_A103)) + geom_histogram(fill="blue", alpha=0.4) + ggtitle("distribution du nombre de passages") # alpha gère la transparence
Comparons cette distribution par type de structure :
p <- ggplot(urgences, aes(x=ID_A103, fill = ID_A4)) + geom_histogram(alpha=0.4) + theme(legend.position="top") # en précisant la couleur par un nom de variable, on demande le tracé d'un histogramme pour chaque modalité
p
Avec la fonction ggplotly de plotly on peut rendre le graphique interactif.
library(plotly)
ggplotly(p)
On va désormais transformer cette variable en variable qualitative pour rassembler les établissements de taille similaire
valeurs_manquantes <- is.na(urgences$ID_A103) # on regarde si la variable admet des valeurs manquantes,
head(valeurs_manquantes)
## [1] FALSE FALSE FALSE FALSE FALSE FALSE
(TRUE %in% valeurs_manquantes) # c'est le cas
## [1] TRUE
sum(valeurs_manquantes) # nb de valeurs manquantes
## [1] 1
dim(urgences)
## [1] 734 73
urgences <- na.omit(urgences) # si on veut retirer toutes les lignes avec au moins une valeur manquante (pas toujours souhaitable)
dim(urgences)
## [1] 725 73
# on crée une nouvelle colonne
urgences$nb_passages_3 <- cut(urgences$ID_A103, breaks = c(0,40,80,max(urgences$ID_A103, na.rm=T))) # max, comme la plupart des fonctions renverra une erreur si appliqué sur un vecteur contenant des NA, on les enlève pour le calcul avec l'option na.rm
table(urgences$nb_passages_3)
##
## (0,40] (40,80] (80,263]
## 187 309 229
On voudrait maintenant regarder le nombre de patients ou les effectifs par tranche horaire, on voit bien qu’il serait plus facile d’avoir une variable tranche horaire que de considérer des colonnes différentes contenant tantôt un A, un B, un C etc… On va récupérer les colonnes correspondant à chaque horaire, les renommer, puis de les empiler en créant une colonne horaire
horaire_8 <- urgences[,c(paste('ID_','A',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2')] # sélection des colonnes contenant un A pour chaque établissement
horaire_12 <- urgences[,c(paste('ID_','B',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2')]
horaire_18 <- urgences[,c(paste('ID_','C',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2')]
horaire_22 <- urgences[,c(paste('ID_','D',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2')]
horaire_8_ <- urgences[,c(paste('ID_','E',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2')]
head(horaire_8)
#on rajoute une colonne d'heure
horaire_8$heure <- '8h'
horaire_12$heure <- '12h'
horaire_18$heure <- '18h'
horaire_22$heure <- '22h'
horaire_8_$heure <- '8h_'
#pour pouvoir empiler les 5 tables il faut qu'elles aient les mêmes noms de colonnes
names(horaire_8) <- c(paste('ID_',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2','heure')
names(horaire_12) <- c(paste('ID_',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2','heure')
names(horaire_18) <- c(paste('ID_',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2','heure')
names(horaire_22) <- c(paste('ID_',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2','heure')
names(horaire_8_) <- c(paste('ID_',c(87:90,92:101), sep =''),'ID_A103','ID_A4','nb_passages_3','ID_A2','heure')
urgences2 <- rbind(horaire_8, horaire_12, horaire_18, horaire_22, horaire_8_) # on les empile avec rbind
dim(urgences)
## [1] 725 74
dim(urgences2)
## [1] 3625 19
head(urgences2)
# on va remplacer les variables 'ID_A103', 'ID_A4' et 'ID_A2' par des noms de variables intelligibles
names(urgences2)[names(urgences2)=='ID_A103'] <- 'nb_passages'
names(urgences2)[names(urgences2)=='ID_A4'] <- 'type_structure'
names(urgences2)[names(urgences2)=='ID_A2'] <- 'index_structure'
On peut regarder le nombre de médecins urgentistes ID_92 en fonction de la tranche horaire et de la taille de la structure avec un graphique à boites.
sapply(urgences2, class) #attention certaines variables numériques ne sont pas considérées comme telles !
## ID_87 ID_88 ID_89 ID_90
## "integer" "integer" "integer" "integer"
## ID_92 ID_93 ID_94 ID_95
## "factor" "factor" "integer" "factor"
## ID_96 ID_97 ID_98 ID_99
## "integer" "factor" "integer" "integer"
## ID_100 ID_101 nb_passages type_structure
## "integer" "integer" "integer" "factor"
## nb_passages_3 index_structure heure
## "factor" "factor" "character"
urgences2[,1:15] <- lapply(urgences2[,1:15],function(v) as.numeric(as.character(v))) # on les convertit
p<-ggplot(data = urgences2, aes(x=heure, y=ID_92, color=nb_passages_3)) + geom_boxplot() + ylab('Effectif de médecin urgentiste')
p
Si on veut plutôt regarder l’évolution du nombre moyen de médecins urgentistes dans des établissements de taille similaire au cours de la journée, il va falloir d’abord calculer cette moyenne par groupe (les groupes étant identifiés par les 3 classes construites précédemment en rendant qualitative la variable des nombres de passages et par l’heure de pointage). Pour cela on va faire appel au package dplyr, vous devez l’installer si ce n’est pas déjà fait.
library(dplyr)
grouped <- group_by(urgences2, heure, nb_passages_3)
grouped_mean <- summarise(grouped, mean=mean(ID_92, na.rm = T))
head(grouped_mean)
# les deux actions peuvent être réalisées dans la foulée avec l'opérateur %>%
grouped_mean <- urgences2 %>% group_by(heure, nb_passages_3) %>% summarise(mean=mean(ID_92, na.rm = TRUE))
head(grouped_mean)
ggplot(data = grouped_mean, aes(x=heure, y=mean, group = nb_passages_3, color = nb_passages_3)) + geom_line() + geom_point()
Remarque : les fonctionnalités du packages dplyr sont bien plus riches, allez jeter un oeil à la cheatsheet data wrangling french fournie dans le dossier de la formation.
Maintenant on va regarder l’évolution de la distribution des effectifs de personnels dans le temps, sans se limiter aux médecins urgentistes. Pour présenter les différents graphiques côte à côte, on va utiliser l’option facet de ggplot2 mais avant cela, il faut retravailler encore la forme des données. Commençons par regarder les distributions d’effectifs de professionnels, donc les variables de 90 à 101, que l’on va rendre plus intelligible dans un premier temps.
urgences_ps <- urgences2[,c('index_structure','nb_passages_3','heure','ID_90','ID_92','ID_93','ID_94','ID_95','ID_96','ID_97','ID_98','ID_99','ID_100','ID_101')]
names(urgences_ps)[4:length(urgences_ps)] <-c('med urgentistes', 'med non urgentistes','internes','médecins int et ext', 'cadres','IDE','AS','ASH','secrétaires','brancardiers') # on rend les colonnes plus intelligibles
library(reshape)
urgences_ps <- melt(urgences_ps, id=c("index_structure","nb_passages_3",'heure'))
head(urgences_ps)
names(urgences_ps)[4:5] <-c('PS','effectif')
head(urgences_ps)
On a donc créé une variable PS (pour professionnel de santé) qui prend plusieurs modalités et on a, pour chaque structure d’urgence identifiée par l’identifiant ID_A2, une ligne par profession avec les effectifs correspondants. Cette représentation rend la représentation graphique quasi immédiate avec ggplot2
urgences_ps$heure = factor(urgences_ps$heure , levels=c("8h", "12h", "18h", "22h", "8h_")) # on transforme la variable en facteur et on ordonne les heures
p <- ggplot(urgences_ps, aes(x=PS, y=effectif, fill=PS)) + geom_bar(stat = "identity") + facet_grid(as.factor(heure) ~ nb_passages_3) + theme( axis.ticks.x=element_blank(),axis.text.x=element_blank())
p
# theme permet de retirer la légende sous l'axe des x, si vous retirez ce bout de code vous verrez la différence
# facet_grid permet de produire plusieurs graphiques en fixant une variable ici le nombre de passage ou l'heure de pointage
Le même exercice peut être réalisé pour la distribution des types de patients, à vous de jouer
Enfin, nous allons représenter la répartition géographique des structures d’urgence en France. Pour cela nous allons utiliser la table sas ej.sas7bdat qui relie les codes FINESS des établissements avec leurs coordonnées géographiques. Cette base est stockée dans un format SAS, on va utiliser la librairie haven pour la lire.
library(haven)
ej <- read_sas(data_file = 'ej.sas7bdat')
head(ej)
Pour tracer sur une carte les structures d’urgence de notre table urgences, on va apparier notre table avec ej, sur la clé finess et ID_A2 respectivement. On va bien faire attention à garder tous les codes finess de la base urgences
urgences_loc <- merge(urgences[,c('ID_A2','ID_A4','ID_A103')],ej,by.x = 'ID_A2',by.y = 'finess', all.x = TRUE) #on garde aussi le type de structure et le nombre de passages
head(urgences_loc)
dim(urgences)
## [1] 725 74
dim(urgences_loc)
## [1] 725 15
dim(ej)
## [1] 55511 13
Le fait d’avoir indiqué l’option all.x = TRUE force à garder toutes les lignes de urgences même quand il n’existe pas de match possible. Ici certains codes finess ne semblent pas avoir de coordonnées géographiques disponibles. On ne va conserver que les structures qui ont des coordonnées grâce à na.omit
urgences_loc_sans_na <- na.omit(urgences_loc)
dim(urgences_loc_sans_na)
## [1] 329 15
Seulement la motié des structures sont géocodées.
On va regarder la répartition des structures d’urgence par type de structure (repéré par la variable ID_A4). La représentation de points sur une carte nécessite de maîtriser le concept de projection : https://medium.com/@_FrancoisM/introduction-%C3%A0-la-manipulation-de-donn%C3%A9es-cartographiques-23b4e38d8f0f.
library(leaflet)
library(rgdal)
urgences_loc_sans_na <- na.omit(urgences_loc)
# Ce bloc permet de convertir les coordonnées en Lambert 93 en WGS 84 plus traditionnel
coordinates(urgences_loc_sans_na) <- c("X", "Y")
proj4string(urgences_loc_sans_na) <- CRS("+init=epsg:2154") # WGS 84
CRS.new <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
urgences_loc_sans_na <- spTransform(urgences_loc_sans_na, CRS.new)
urgences_loc_sans_na$lon <- data.frame(coordinates(urgences_loc_sans_na))$X
urgences_loc_sans_na$lat <- data.frame(coordinates(urgences_loc_sans_na))$Y
factpal <- colorFactor(topo.colors(3), urgences_loc_sans_na$ID_A4) # palette de 3 couleurs
m <- leaflet(urgences_loc_sans_na) %>% addTiles() %>%
addCircles(lng = ~lon, lat = ~lat, weight = 1,
radius = ~ID_A103*100,
popup = ~paste(rs, ", nb de passages : ", ID_A103, ", type de structure :", ID_A4),
color = ~factpal(ID_A4), fillOpacity = 0.5) %>%
addLegend("bottomleft", title = "Type de structure d'urgences", pal = factpal, values = ~ID_A4, opacity = 0.7)
m